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Abstract.  A  procedure  is  developed  that,  in  two  iterations,  solves  the 

hyperbolic  Kepler's  equation  in  a  ver>  efficient  manner,  and  to  an 

-20 

accuracy  that  proves  to  be  always  better  than  10  (relative  truncation 
error).  Earlier  work  on  the  elliptic  equation  has  been  extended  by  the 
development  of  a  new  procedure  that  solves  to  a  maximum  relative  error 
of  10~14. 

1.  Introduction 

In  an  earlier  paper  (Odell  and  Gooding,  1986),  we  considered  the 
classical  equation  of  Kepler  for  elliptic  orbits,  and  recommended  two 
particular  procedures  for  its  solution.  This  paper  was  a  greatly 
shortened  version  of  the  authors'  monograph  (Gooding  and  Odell,  1985), 
and  all  references  to  our  'previous  work'  in  the  present  paper  should  be 
taken  to  cite  the  two  references  just  given.  To  complement  the  previous 
work,  we  have  now  studied  the  hyperbolic  equation,  and  the  present  paper 
is  a  shortened  version  of  the  recent  RAE  report  by  Gooding  (1987). 

Though  there  is  a  very  extensive  literature  on  the  elliptic  equation, 
several  papers  having  appeared  since  our  previous  work,  little  has  been 
published  on  the  hyperbolic  equation.  Burkardt  and  Danby  (1983)  con¬ 
sider  it  briefly,  before  passing  to  a  generalized  equation  expressed 
in  terms  of  universal  variables,  whilst  other  authors,  such  as  Bergam 
and  Prussing  (1982),  only  consider  hyperbolic  orbits  in  the  context  of 
the  universal  (generalized)  equation. 

The  use  of  universal  variables  is  aesthetically  very  attractive 
and  we  discussed  the  universal  equation  briefly  in  our  previous  work, 
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referring  in  particular  to  the  recent  paper  by  Shepperd  (1985).  Since 
then,  one  of  us  (Gooding,  1989)  has  been  studying  the  use  of  universal 
elements  in  conversion  algorithms  (to  and  from  position  and  velocity), 
but  a  conclusion  of  that  study  is  that  the  best  way  to  proceed,  inside 
the  relevant  algorithm,  is  still  to  solve  the  elliptic  equation,  the 
hyperbolic  equation,  or  Barker's  equation,  as  appropriate. 

The  hyperoolic  Kepler’s  equation  is 

e  sinh  H  -  H  =  ,  (1) 

where  e  (eccentricity)  and  (hyperbolic  mean  anomaly)  are  assumed 

known,  and  H  (hyperbolic  eccentric  anomaly,  often  denoted  by  F  )  is 
to  be  determined.  Since  the  range  of  e  is  from  unity  (rectilinear 
hyperbolic  orbit)  to  infinity  (uniform  linear  motion,  for  zero  attrac¬ 
tive  force),  there  is  some  advantage  in  rewriting  equation  (1)  as 

sinh  H  -  gK  »  L  ,  (2) 

where 

g  «  1/e  ,  L  »  M^/e  ; 

we  also  define 

8X  "  1  -  g  • 

Later  we  shall  find  further  advantage  in  reformulating  the  equation  so 
that  sinh  H  ,  rather  than  H  ,  is  the  quantity  to  be  determined. 

Our  approach  to  the  iterative  solution  of  equation  (1),  or  equiv¬ 
alently  (2),  has  been  to  adhere  as  far  as  possible  to  the  philosophy 
underlying  the  final  procedure  that  we  previously  recommended.  In  par¬ 
ticular,  we  have  sought  a  combination  of  starting  formula  and  iteration 
process  accurate  enough  for  the  resulting  procedure  to  have  always  met  a 


5 


given  accuracy  goal  after  a  fixed  small  number  (preferably  two)  of 
iterations;  assuming  a  smoothly  continuous  starter  over  the  (e ,  M^)- 
data  space,  the  output  (H  or  sinh  H)  would  then  also  be  smoothly  con¬ 
tinuous.  ('Smooth  continuity'  is  meant  to  be  synonymous  with  the  term 
'smooth  portability'  discussed  in  the  previous  work.)  The  accuracy  goal 

was  expressed  as  a  ceiling  for  the  relative  error  in  the  output  quan- 
-13 

tity,  the  ceiling  of  10  being  selected;  as  in  the  previous  work,  this 
was  related  to  the  computer  used  for  most  of  the  work,  but  the  goal  has 

not  been  the  same,  since  for  the  elliptic  equation  the  value  of 
-13 

10  (rad)  had  been  an  absolute-error  ceiling. 

Section  2  describes  two  procedures  that  solve  for  H  ,  and 
Section  3  describes  two  procedures  that  solve  for  S  (  ”  sinh  H).  As 
the  previously-developed  iterator  (described  as  'HN'  and  providing 
quartic  convergence)  was  incorporated  in  all  four  procedures,  Sections  2 
and  3  mainly  provide  descriptions  of  the  starters.  Section  A  is  devoted 
to  computer  implementation  and  results  for  the  last  (and  best)  of  the 
four  procedures,  and  Section  5  describes  a  new  procedure  developed  for 
the  elliptic  equation,  such  that  the  relative-error  goal  used  for  the 
hyperbolic  equation  can  be  met.  Also  included  in  Section  5  is  a 
discussion  of  an  alternative  iteration  method,  due  to  Laguerre  and 
recently  recommended  by  Conway  (1986);  though  it  only  provides  cubic 
convergence,  it  is  much  more  robust  than  our  standard  iterator,  so  it 

comes  into  its  own  when  a  good  starter  is  not  available. 

-13 

The  recommended  hyperbolic  procedure  meets  the  10  goal  with  the 
greatest  of  ease,  since  it  provides  20-decimal-digit  accuracy  in  all 
cases.  The  new  elliptic  procedure  does  not  do  so  well,  the  accuracy 
provided  being  1A  digits. 


6 


2.  Procedures  solving  for  H 


Our  initial  approach  was  to  follow  EKEPT.2  (the  final  procedure  developed 
in  the  previous  work)  as  closely  as  possible.  The  starting  formula  for 
EKEPL2  was  based  on  an  interpolation  between  starters  suitable  for 
e  ”  0  and  e  =  1  ,  the  latter  (Eq^)  being  a  composite  starter  smoothly 
patching  together  a  cube-root  expression  and  a  bilinear  expression. 

Thus  to  solve  (2),  we  required  a  special  starter,  ,  for  g  «  1  , 

after  which  the  overall  starting  formula  would  be  given  by 

H0  -  8  H01  +  8l  slnh  1  L  .  O) 

since  sinh  ‘L  is  the  limiting  solution  of  (2)  for  g  ■*  0  .  We  assume, 
for  convenience,  that  is  positive  and  that,  as  in  the  previous 

work,  HQ1  has  a  pair  of  components  that  are  to  fit  smoothly  together  at 
Mh  »  1/6  .  Then  a  suitable  definition  is 


I ( 6Mh) 1/3  if  0  <  <  1/6  (4a) 

In  2(1^  +  1/3)  +  1  if  1/6  .  (4b) 


When  the  overall  starter  was  coupled  to  the  standard  iteration 
process,  however,  the  resulting  procedure  gave  somewhat  disappointing 
results.  In  particular,  the  accuracy  after  two  iterations  was  much 
worse  than  for  EKEPL2.  As  shown  by  the  detailed  analysis  of  the 
(hyperbolic)  procedure  by  Gooding  (1987)  ^  in  fact,  the  guaranteed 
accuracy  is  only  seven  decimal  digits  (relative  truncation  error),  so  it 
was  necessary  to  extend  the  procedure  to  a  third  iteration.  This  led  to 
an  accuracy  of  at  least  26  digits  guaranteed,  compatible  with  expec¬ 
tation  on  the  basis  of  quartic  convergence,  but  the  need  for  a  third 
Iteration  detracted  from  the  efficiency  of  the  procedure.  The  procedure 
was  in  any  case  less  efficient  than  EKEPL2,  since  equation  (4b)  cannot 


be  evaluated  as  cheaply  as  a  bilinear  expression.  But  it  was  a  third 
defect  of  the  procedure  that  was  of  most  interest,  and  this  is  described 
in  the  next  paragraph;  the  defect  is  also  suffered  by  EKEPL2,  which  is 
why  further  attention  was  in  due  course  given  to  the  elliptic  equation 
as  well  as  the  hyperbolic. 

For  small  values  of  ,  equation  (2)  gives 

H3  +  6glH  -  6L  -  0  ,  (5) 

and  this  is  the  basis  for  (4a),  which  is  clearly  an  excellent  starter 
when  g^  *  0  .  For  g^  >  0  ,  however,  the  neglect  of  the  second  term  of 
(5)  can  lead  to  a  serious  overestimate  of  H  ,  which  the  interpolation 
by  (3)  does  little  to  remedy  when  g^  is  only  slightly  greater  than 
zero  (e  only  slightly  greater  than  unity).  The  overestimation  is  never 
so  gross  as  to  affect  convergence  adversely,  so  long  as  this  statement 
is  interpreted  in  terms  of  truncation  error  (and  this  applies  to  rela¬ 
tive  error,  not  just  absolute  error),  but  the  effect  on  rounding  error 
is  another  matter  entirely,  as  shown  in  the  detailed  analysis.  Even 
here  there  is  no  problem  with  absolute  error,  because  the  iterator  works 
so  well,  but  relative  rounding  error  can  only  be  reduced  at  a  rate  (per 
iteration)  deLerminad  by  the  computer  word-length. 

As  already  remarked,  the  defect  of  large  relative  rounding  error 
applies  to  EKEPL2,  from  the  previous  work,  as  well  as  to  the  hyperbolic 
procedure  being  described.  It  was  not  observed  before,  however,  because 
the  previous  analysis  was  conducted  almost  entirely  in  terms  of  absoluze 
error,  this  being  possible  because  of  the  periodic  nature  of  E  -  M  for 
the  elliptic  equation;  thus  absolute  error  is  bounded  over  the  full 
(infinite)  range  of  M  .  In  consequence,  we  did  not  look  at  very  small 
non-zero  values  of  M  ;  had  we  done  so,  the  problem  would  have  come  to 
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light  In  the  plotting  of  rounding  error  in  Fig  6  of  Odell  and  Gooding 
(1986). 


With  a  view  to  eliminating  the  defect  just  described,  and  in  the 
hope  of  devising  a  procedure  for  which  two  iterations  would  always 
suffice,  we  abandoned  our  initial  approach  to  the  solution  of  (2),  in 
favour  of  one  based  on  the  approximation  (5).  The  new  approach,  as  in 
our  previous  consideration  of  the  approximating  cubic  for  the  elliptic 
equation,  follows  the  notation  of  Ng  (1979). 

We  define  Hqq  as  the  starter,  which  is  a  good  approximation  for 
small  and  arbitrary  g  ,  given  by  the  (unique)  solution  of  the 

equation  (cf  (5)) 

Hoo  +  3«  Hoo  -  2r  -  0  *  (6> 

where  q  *  2g^  and  r  -  3L  .  As  before,  we  note  that  (6)  is  the 
classical  cubic  equation,  for  which  we  can  do  better  than  Ng  by 
expressing  the  solution  with  a  single  cube  root,  as 

Hon  -  s  -  q/s  »  C7) 

where  (assuming  ,  and  hence  r  ,  to  be  non-negative) 

s  ■  |V +  +  r] 

This  expression  for  Hqq  is  still  not  optima’,  however,  as  unnecessary 

2 

rounding  error  will  be  experienced  when  s  *  q  ;  some  algebraic 
manipulation  gives  us  the  optimal  solution  as 
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basis  for  an  overall  starter.  Thus  Is  defined,  by  use  of  weights 

and  1,  as 

Ho  ’  (Mh  slnh_lL  +  Hoo)/(‘\  +  •  (9) 

On  coupling  equation  (9)  to  Che  usual  iterator,  we  have  a  pro¬ 
cedure  that  Is  free  of  the  deficiencies  of  the  first  procedure.  In 
particular,  there  is  no  rounding-error  problem,  and  the  maximum  relative 
error  (due  to  truncation)  after  two  iterations  is  only  7  x  10 

3 .  Procedures  solving  for  S  (  •  sinh  H) 

The  development  of  the  second  procedure  of  Section  2  might  have  ended 
our  study  of  the  hyperbolic  equation,  but  an  alternative  approach  was 
thought  worth  looking  at.  In  this,  equation  (1)  is  not  just  rewritten 
as  (2),  but  more  significantly  rewritten  as 

S  -  g  sinh  1 S  ~  L  ,  (10) 

where 

S  -  sinh  H  .  (11) 

Then  S  ,  instead  of  H  ,  is  the  quantity  to  be  solved  for. 

There  are  at  least  two  advantages  in  solving  equation  (10),  as 
opposed  to  (1)  or  (2).  First,  it  is  more  efficient,  since  the  only  use 
to  be  made  of  H  is  normally  via  sinh  H  and  cosh  H  ;  moreover,  there 
is  no  loss  of  accuracy  in  deriving  cosh  H  as  /(I  +  S2 )  ,  in  complete 
contrast  to  the  loss  that  can  occur  if  cos  E  is  derived  from  sin  E  . 
Secondly,  there  is  a  gain  in  precision  with  the  use  of  S  ,  rather  than 
H  ,  for  large  values  of  J  S  J  ;  thus  the  evaluation  of  the  composite 
slnh(sinh  function,  like  the  composite  exp(ln)  ,  cannot  be  relied 


upon  to  be  the  identity  function  when  the  argument  is  large.  It  is  to 
be  noted,  also,  how  great  is  the  resemblance  between  the  reformulated 
equation  and  the  standard  elliptic  equation,  since  g  has  the  same 

range  as  e  for  an  ellipse,  whilst  sinh  ,  like  sin  x  ,  expands  as 

3  ,5. 

x  -  x  /6  +  0[x  ]  .  As  a  final  introductory  point,  we  remark  thac  the 
change  of  variable  from  H  to  S  is  bound  to  affect  the  operation  of 
the  iteration  process;  thus  if  we  chose  to  solve  both  equations  (1)  and 
(10),  with  starters  (H  and  S  )  related  by  (11),  this  relation  would 
not  be  preserved  during  iteration. 

The  first  approach  to  a  starter  for  equation  (10)  was  based  on  the 
direct  modification  of  the  successful  solution  of  equation  (2)  via  the 
cubic  equation  (6).  Thus  we  solve  (6)  for  Sqq  ,  rather  than  Hqq  ,  if 
we  replace  Che  definitions  of  q  and  r  by 

q  -  2gl/g  and  r  -  3l/g  .  (12) 


In  view  of  the  possibly-zero  denominators  in  (12),  however,  it  is  better 
if  we  write  the  cubic  equation  in  Che  homogeneous  form 


a  S„  -  +  3b  -  2c 


"00 


00 


(13) 


where  a  ■  g  ,  b  "  2g^  (  »  gq)  and  c  »  3L  (  *=  gr)  . 

From  Sq0  ,  the  starter  appropriate  to  small  L  (still,  for  con¬ 
venience,  assumed  positive),  we  require  the  extension  to  an  overall 
starter,  in  analogy  with  (9).  The  most  obvious  analogy  seemed  to  be 
with  Sq  given  by 

s0  -  (l2+soo)/(l+  °  • 

but  this  gave  disappointing  results.  Afte”  some  experimentation,  it 
found  that  the  starter  defined  by 

SQ  -  SQ0  +  LZ/(L  +  5)  (U) 


was 


!  ] 


worked  much  better,  the  constant,  5,  having  been  determined 
empirically.  The  results  remained  worse,  rather  than  better,  than  chose 
derived  using  the  second  procedure  for  H  ,  however,  so  it  was  decided 
to  investigate  an  entirely  different  type  of  starter. 

One  of  the  possible  starters  (the  eleventh)  discussed  in  our 
previous  work  was  a  very  complicated  one,  constructed  on  Che  rationale 
that  it  should  have  the  right  form  in  the  region  of  awkward  convergence, 
whilst  at  the  same  time  matching  the  formal  e-series  expansion  (of  the 
solution  to  the  elliptic  equation)  to  terms  of  order  e3 .  By  'having 
the  right  form'  etc,  we  just  meant  that,  for  e  -  1  and  small  M  ,  the 
starter  should  behave  like  V 6M  ,  since  the  more  stringe  it  constraint 
(Imposed  by  the  cubic-equation  approximation,  cf(5),  for  e  *  1)  would 
then  be  met  by  the  series  matching.  It  was  decided  to  apply  this 
rationale  to  (the  hyperbolic)  equation  (10),  but  to  develop  the  desired 
new  starter  in  a  straightforward  and  systematic  manner  by  basing  it  on 
Lagrange's  expansion  theorem  (Whittaker  and  Watson,  1940). 

The  objective  may  be  summarized  as  being  an  expression  of  the  form 

SQ  -  L  +  g  B_1/3  sinh-1L  ,  (15) 

2 

where  B  *  1  +  terms  in  g  ,  g  ,  etc  as  required.  Clearly 

B  -  F(S)/F(L)  ,  (16) 

,  -3 

where  F  is  the  function  (sinh  )  (a  somewhat  unfortunate  notation  as 
'he  tv  >  negative  signs  have  different  meanings).  We  will  use  the 
xpanslon  theorem  to  express  F(S)  as  a  function  of  L  (with  parameter 
g  ),  and  to  help  in  this  we  define  l  *  sinh  3L  .  Then  the  theorem 


gives 


F(S) 


(H) 


-  F(L)  +  ^  |y  D(n_1)(inD  F(L))  , 

n-1 

where  D  is  the  operator  d/dL  . 

On  performing  the  necessary  differentiations  and  dividing  by 
F(L)  ,  we  find  that  (17)  leads  to 

B  -  (1  -  gA]3  +  ^g2lA3(3L(l  -  gA)  +  gAA2(l  -  2L2))  +  o(g4)  ,  (18) 

2 

where  A  »  sech  l  .  But  B  has  to  behave  like  l  /6  ,  when  g  *  1  and 
i  is  small;  so  we  give  up  the  formally-accurate  g3  term  in  (18), 
writing  instead 

B  -  (1  -  gA)3  +  g2f.LA3(9  -  8g)/6  +  0(g3)  .  (19) 

Substitution  of  (19)  in  (15)  gives  the  required  starter.  Coupled 
with  the  usual  Iteration  process,  it  works  so  well  (as  we  will  see  in 
the  next  section)  that  the  resulting  procedure  constitutes  our 
unqualified  recommendation  for  use  in  solving  the  hyperbolic  Kepler's 
equation. 

4.  Implementation  and  Results  for  the  Recommended  Procedure 

The  Fortran-77  function  SHKEPL  has  been  written  to  implement  the 
recommended  solution  procedure.  Its  arguments  are  L  and  g^  ,  and  it 
is  listed  in  Appendix  A.  The  second  argument  is  g^  ,  rather  than  g  , 
to  minimize  rounding  error  in  the  vicinity  of  (g,L)  ”  (1,0),  when 
S  -  g  sinh  has  to  be  computed  by  the  special  subordinate  function 
SHMKEP,  listed  in  Appendix  B;  SHMKEP  operates  in  the  tame  way  as  the 
function  (EMKEPL)  we  gave  before,  except  that  it  is  rather  more 
efficient  since  S  -  g  sinh  is  computed  as  g^S  +  g(S  -  sinh  *S) 
rather  than  (S  -  sinh  ^S)  +  g^  sinh  ''S  •  Two  other  (non-standard) 
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subordinate  functions  are  used  by  SHKEPL:  DCUBRT  for  obtaining  a  cube 
root;  and  DASINH  for  the  inverse  sinh  • 

The  computer  used  for  testing  SHKEPL  has  mainly  been  a  PRIME  750, 
as  in  the  previous  study,  but  a  Cray  IS  was  used  for  the  most  accurate 
results.  Though  the  input  for  SHKEPL  is  actually  L  and  g^  ,  the  test 
data  consisted  of  a  wide  range  of  values  for  H  and  e  ,  H  being  the 
source  for  the  reference  value  of  S  ,  from  which  the  nominally  input  L 
was  obtained  via  equation  (10).  From  the  value,  S^  ,  extracted  from 
SHKEPL  after  i  iterations,  we  have  the  relative  error,  ,  given  by 

cr1  -  (S1  -  S)/S  .  (20) 

-30  4 

For  each  e  ,  the  range  of  H  extended  from  10  to  10  ,  at  a 
fixed  geometric  interval,  the  maximum  value  of  |  J  obtained  being  a 
measure  of  the  accuracy  of  SHKEPL  with  i  iterations.  (For  normal  use, 
and  as  listed  in  Appendix  A,  1  is  fixed  at  2,  but  the  values  0  and  1 
were  also  of  Interest.)  A  very  wide  range  of  e  was  tested,  extending 
from  unity  to  10^®. 

Figure  l  provides  plots  of  log  j  Pmax  |  against  log  e  ,  for 
1  *  0  (starter),  1  and  2  (standard  SHKEPL).  The  range  of  e  ,  here, 
extends  only  from  1  to  10,  as  the  accuracy  is  so  good  for  higher  values 
of  e  .  The  accuracy  is,  in  fact,  striking,  since  it  is  evident  that 
SHKEPL  gives  at  least  20  significant  figures  correct  (for  S)  in  all 
circumstances.  If  the  function  is  restricted  to  a  single  iteration,  the 
number  of  significant  figures  reduces  to  five  -  this  is  compatible  with 
the  quartic  nature  of  the  iterator. 

5.  The  Elliptic  Equation  revisited 

In  our  previous  study  we  established  the  merits  of  two  particular 
procedures  for  solving  the  standard  Kepler  equation 


E  -  e  sin  E 


M 


(21) 
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The  second  (and  more  accurate)  of  these  procedures  (EKEPL2)  combined  a 
number  of  desirable  features,  as  described  in  our  papers;  in  particular, 
the  incorporation  of  a  bilinear  formula  in  the  starter  made  it  extremely 
efficient.  We  failed  to  observe,  however,  that  EKEPL2  suffers  from  the 
relative-rounding-error  defect  remarked  on  for  the  first  procedure  of 
Section  2.  This  may  be  regarded  as  a  very  minor  weakness,  which  it  does 
not  seem  possible  to  remedy  efficiently;  for  completeness ,  however, 
we  describe  a  procedure,  modelled  on  the  second  procedure  of  Section  2, 
that  is  free  of  this  defect  and  retains  all  the  merits  of  EKEPL2  other 
than  efficiency.  We  also  take  the  opportunity  to  make  some  comments  on 
the  application  of  Laguerre  iteration  to  Kepler's  equation,  in  the 
light  of  the  recent  paper  by  Conway  (1986). 

We  have  seen,  in  the  context  of  the  hyperbolic  equation,  that  the 
relative-rounding-error  defect  could  be  eliminated  by  use  of  a  cubic 
approximation,  for  small  H  ,  valid  for  all  e  and  not  just  e  «  1  . 

We  adopt  the  same  policy  for  the  elliptic  equation,  requiring  to  solve 
the  equation,  of  (6) , 

E0 1  +  ^00  -  2-  -  0  ,  (22) 

given  in  our  previous  work,  where  now 

q  «  2e^/e  and  r  “  3M/e  ,  with  =  1  -  e 

Again  we  prefer  to  solve  the  homogeneous  equation,  parallel  to  (13);  the 
function,  DCBSOL,  for  doing  this  has  (like  DCUBRT  and  DASINH)  been  listed 
by  Gooding  ( 1 987)  . 
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The  solution  of  (22)  for  E^q  provides  a  suitable  starter  for 
small  M  (>  0)  and  any  e  (0  <  e  <  1)  ,  so  for  arbitrary  M  over  the 
range  [0,ir]  we  use  the  formula 


E 


0 


M2  +  (r  -  M)Eoq 

TT 


(23) 


Coupled  to  our  usual  iterator,  this  led  to  a  new  procedure,  EKEPL3,  of 
arguments  (as  for  EKEPL1  and  EKEPL2)  M  and  e  .  To  minimize  rounding 
error  in  the  awkward  region  where  (e ,M)  •  (1,0),  the  second  argument  was 
then  changed  to  e^  and  the  resulting  procedure  given  the  name  EKEPL. 
The  Fortran-77  function  EKEPL  is  listed  in  Appendix  C;  it  uses  the 
subordinate  function  EMKEP  in  the  awkward  region,  but  this  function  is 
so  like  the  original  EMKEPL  that  it  is  not  listed  here. 

Figure  2  provides,  for  the  solution  procedure  EKEPL,  plots  of 


Where  «  (Ei  -  E)/E  ,  (24) 

for  a  (decreasing)  range  of  e  from  unity  to  10  0-275  (»  0.253).  To 
provide  a  comparison  with  the  old  procedure  (EKEPL2),  the  corresponding 
curve  for  that  procedure  is  also  given,  for  the  definitive  i  -  2  .  It 
is  seen  that  the  new  procedure  does  better,  except  when  e  >  0.95 
(log  e  >  -0.023)  ,  and  there  is  an  apparent  paradox  here,  since  it  was 
for  precisely  the  high  values  of  e  (approaching  unity)  that  EKEPL2  was 
defective.  The  paradox  is  easily  resolved,  however,  since  the  defect 
relates  to  rounding  error,  whereas  Figure  2  only  covers  truncation 
error:  for  e  *  l  ,  both  procedures  give  very  small  truncation  error 

when  M  is  small  (ie  in  the  'awkward  region'!)  and  the  curves  of 
Figure  2  are  determined  by  'macroscopic'  values  of  H  ;  it  is  a  tribute 


to  the  efficacy  of  the  bilinear  function,  used  in  EKEPL2,  that  under 
these  circumstances  EKEPL2  does  so  much  better  than  the  new  EKEPL. 

Figure  2  indicates  that  the  accuracy  for  EKEPL  is  never  worse  than 
about  14  significant  figures.  Though  good,  this  is  inferior  to  the 
accuracy  of  SHKEPL.  Since  SHKEPL  uses  a  starter  based  on  Lagrange's 
expansion  theorem,  it  might  be  thought  that,  as  for  the  reformulated 
hyperbolic  equation,  we  would  do  much  better  via  this  theorem  than  via 
the  cubic-equation-based  starter.  This  is  not  so,  however;  results  for 
a  procedure  based  on  the  expansion  theorem,  derived  in  exactly  the  same 
way  as  for  the  hyperbolic  equation,  are  actually  worse  than  those  given 
by  EKEPL.  Despite  this,  we  give  the  equations  corresponding  to  (15)  and 
(19),  for  completeness;  they  are 

Eq  »  M  +  eB_1/3  sin  M 

and 

B  *  (l  -  e  cos  M)3  +  e2  sin  M  (9  -  8e)/6  +  0(e  ) 

We  now  discuss  the  merits  of  replacing  our  quartic  generalized- 
Newton  iteration  formula  by  one  of  the  Laguerre  formulae,  which,  though 
they  only  give  cubic  convergence,  are  shown  by  Conway  (1986)  to  be  much 
more  robust  than  the  Newton  formulae. 

The  general  Laguerre  formula,  for  iterating  to  a  root  of  the 
equation  f(x)  “  0  ,  may  be  written 

6  -  -  - — -  ,  (25) 

f'jl  +  /[(N  -  l)2  -  N(N  -  l)ff"/f'2]} 

where  5  »  x  -  Xj^  and  f  is  shorthand  for  f(j<i_j),  etc;  N  may  be 
identified  as  the  degree  of  the  polynomial  that  is  matched  to  f  .  Thus 
N  is  normally  taken  to  be  a  small  integer  (>  2  ,  since  for  N  ”  1 
formula  reduces  to  the  Newton-Raphson  formula,  of  only  quadratic 


the 


) 


convergence),  but  it  does  not  have  to  be  restricted  in  this  way  -  the 
formula  is  sound  with  N  large,  infinite,  or  even  non-integral.  Conway 
expresses  the  denominator  of  (25)  as  f'  ±  /[(n  -  l)2f'2  -  N(N  -  1 ) f f ” ] 
this  provides  a  result  in  the  (unrealistic)  circumstance  that  f'  ■  0  , 
but  at  the  expense  of  the  sign  ambiguity  (resolved,  when  f  *  0  ,  by 
giving  the  square  root  the  same  sign  as  f'  )• 

There  would  be  no  advantage  in  substituting  one  of  the  Laguerre 
formulae  for  our  standard  iterator  in  either  of  the  procedures  EKEPL2  or 
EKEPL,  or  in  any  of  the  procedures  developed  for  the  hyperbolic 
equation,  since  the  starters  in  all  these  procedures  are  good  enough  for 
the  full  (quartic)  power  of  the  existing  iterator  to  operate.  For 
EKEPL1  it  is  another  matter,  however,  since  the  possibly  slow 
convergence  here  (infinitely  slow  in  the  limit)  is  completely  overcome 
with  one  of  the  Laguerre  formulae.  The  advantage  of  iteration  by  (25) 
is  even  more  marked  with  the  starter  Eg  •  M  ,  for  which  the  Newton 
formulae  can  lead  to  divergence. 

Conway  demonstrates  the  success  of  the  Laguerre  iterator  in 
solving  Kepler's  equation,  taking  the  arbitrary  value  of  N  =  5  in 
equation  (25).  This  value  is  essentially  a  compromise,  since  in 
straightforward  cases  convergence  is  usually  more  rapid  with  N 
infinite.  In  the  awkward  region  (where  the  Newton  formulae  fail  or  are 
slow),  however,  much  the  best  results  are  given  with  N  •  3  ,  and  the 
reason  for  this  has  been  given  by  Gooding  (1987).  The  essence  of  the 
matter  is  the  following:  for  a  function  approximating  to  the  polynomial 
xn  -  A  (  >  0),  and  an  iteration  step  from  a  gross  overestimate  of  the 
n'th  root,  the  only  value  of  N  ,  in  (25),  that  makes  the  step  a 
quadratic  (as  opposed  to  linear)  one  is  N  -  n  ;  for  Kepler's  equation, 
the  appropriate  value  of  n  is  3,  so  that  N  -  3  is  the  natural 
choice.  What  happens  in  practice,  using  N  ■  3  and  a  poor  initial 


value  (and  assuming  an  awkward-region  problem  still),  is  that  iteration 
steps  are  initially  quadratic;  the  process  then  slows  down  for  a  step  or 
two,  but  finally  homes  in  with  the  full  rapidity  of  its  nominally  cubic 
nature.  When  the  starter  is  Eq  -  M  ,  which  underestimates  E  , 
convergence  seems  to  be  marginally  better  if  the  first  step  is  taken 
with  N  infinite,  before  switching  (because  E^  then  overestimates 
E,  M  being  small)  to  N  ■  3  .  This  is  of  particular  significance, 
since  the  first  step  then  has  the  simple  analytical  outcome 


M  + 


e  sin  M 


/(l  -  2e  cos  M  +  e  ) 


(26) 


The  right-hand  side  of  (26)  is  simply  the  starter  for  EKEPL1;  it  is  an 
attractive  formula,  and  was  first  used  by  Brown  (1931). 

The  conclusions  in  regard  to  the  Laguerre  formulae  are  as  follows. 
When  a  good  starter  is  available,  they  give  no  advantage.  When  the 
starter  is  poor,  on  the  other  hand,  there  is  an  immense  advantage  to  be 
had,  especially  with  N  »  3  .  The  impact  on  our  previous  work  is  that 
the  procedure  EKEPL1  is  greatly  improved  if  the  Halley  iterator  is 
replaced  by  the  Laguerre  iterator  with  N  «  3  . 


6.  Conclusion 


We  have  sought  to  complement  our  previous  work,  on  Kepler's  (elliptic) 
equation,  by  applying  the  same  philosophy  to  the  hyperbolic  equation. 

A  number  of  solution  procedures  were  developed  during  the  new  work,  but 
one  of  these  seems  superior  in  all  respects  and  we  unreservedly 
recommend  it.  It  is  based  on  a  reformulation  of  the  equation,  such  that 
sinh  H  rather  than  H  (hyperbolic  mean  anomaly)  is  determined 
directly;  implemented  in  Fortran-77  under  the  name  SHKEPL,  it  is  listed 
in  Appendix  A. 
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The  procedure  SHKEPL  operates  with  a  starting  formula  derived  by 
use  of  Lagrange's  expansion  theorem  and  with  the  quartic  iteration 
process  developed  in  our  previous  work.  Its  accuracy  is  such  that,  in 
the  absence  of  rounding  error,  two  iterations  (built-in)  should  give  20 
decimal  digits  correct  in  all  cases.  The  effects  of  rounding  error  have 
been  held  to  a  minimum,  such  that  at  most  one  (decimal)  digit  of 
precision  should  be  lost  when  the  computer's  word-length  does  not  exceed 
20  digits. 

During  the  study,  it  was  recognized  that  a  particular  relative¬ 
rounding-error  defect,  identified  in  the  first  procedure  developed, 
would  also  apply  to  the  more  accurate  of  the  two  solution  procedures 
recommended  for  the  elliptic  equation.  Though  it  is  scarcely  con¬ 
ceivable  that  the  defect  could  be  of  consequence  in  practice,  we  have 
developed  an  alternative  procedure,  of  comparable  accuracy  but  not  so 
efficient;  given  the  name  EKEPL,  it  is  listed  in  Appendix  C.  The 
comments  of  the  preceding  paragraph  (on  SHKEPL)  apply  to  EKEPL,  except 
that  the  nominal  worst-case  accuracy  (after  two  iterations)  is  only 
14  digits. 

For  the  elliptic  equation,  which  has  to  be  solved  so  much  more 
often  than  the  hyperbolic  equation,  we  recommend  all  three  solution 
procedures  (the  two  old  ones,  EKEPL1  and  EKEPL2,  and  the  new  one,  EKEPL) 
as  being  appropriate  in  different  circumstances.  A  significant 
Improvement  in  the  least  sophisticated  of  them  (EKEPL1)  can  be  made, 
however,  if  the  Halley  iterator  is  replaced  by  the  Laguerre  iterator 
with  N  =  3  . 

We  have  not  produced  a  'universal  procedure',  for  the  universally 
formulated  Kepler's  equation,  because,  in  spite  of  the  mathematical 
elegance  of  such  a  formulation,  we  regard  it  as  of  little  practical 
value.  For  numerical  work  that  is  accurate  and  efficient,  it  will  always 
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be  necessary  Co  solve  different  equations  for  the  ellipse,  parabola  and 
hyperbola.  However,  this  in  no  way  hinders  the  development  of  outwardly 
universal  algorithms  for  conversion  between  position  and  velocity,  on 
the  one  hand,  and  a  universal  set  of  orbital  elements,  on  the  other; 
such  conversion  algorithms  constitute  the  subject  matter  for  the  com¬ 
panion  paper  (Gooding,  1989). 

In  the  context  of  universal  computation,  we  wish  to  end  with  a 
remark  on  the  recent  response  by  Danby  (1987)  to  our  previous  comment 
(Odell  and  Gooding,  1986)  on  the  disadvantage  of  using  the  Stumpff 
function,  S(x)  ,  when  x  corresponds  to  a  multirevolution  angle  in 
an  elliptic  orbit.  As  he  observes,  the  computation  must  then  be  based 
on  recurrence  formulae,  after  the  angle  has  been  reduced  (by  factors 
of  4)  to  a  suitable  magnitude.  It  is  unfortunate  that  three  of  his 
four  recurrence  formulae,  (16),  are  incorrectly  stated,  the  correct 
versions  being:- 

2 

c0(4x)  =  2  [c o (x) ]  -  1  ,  c i (4x)  »  c0Cx)ci(x)  , 

2 

c2(4x)  =  i  [c  i  (x)  I  ,  c3(4x)  *  {  (c3(x)  +  ci(x)c2(x)J 

Considering  just  the  first  of  these  relations,  we  see  that  when¬ 
ever  co(x)  is  close  to  unity,  in  particular  while  |x|  is  small,  the 
rounding  error  will  be  roughly  quadrupled  in  each  application  of  the 
formula.  The  build-up  will  not  be  as  rapid  as  this  all  the  way  from 
x  »  10  1  to  10^  (to  follow  Danby's  example),  but  the  overall  effect 

will  still  be  the  loss  of  more  than  three  decimal  digits  in  Cq(IO^)  . 

3 

If  cos  E  is  evaluated,  on  the  other  hand,  with  E  =  10  rad  (because 
2 

x  »  E  )  ,  the  intrinsic  loss  (from  storage  of  E  itsalf)  is  at  least 
a  digit  less.  Danby's  technique  is  certainly  viable  (.a  much  more  rapid 
build-up  of  rounding  error  might  have  been  expected  intuitively),  but 
how  much  simpler  and  more  efficient  to  recognize  an  elliptic  orbit  as 
such,  and  apply  old-fashioned  range-reduction! 
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Appendix  A 

THE  SHKEPL  PROCEDURE 

DOUBLE  PRECISION  FUNCTION  SHKEPL  (EL,  Gl) 

EQUATION  EL  =  SHKEPL  +  (Gl  -  1) *DASINH (SHKEPL) , 

WITH  Gl  IN  RANGE  0  TO  1  INCLUSIVE,  SOLVED  ACCURATELY. 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

PARAMETER  (SW=0.5DO,  AHALF=0 . 5D0 ,  ASIXTH=AHALF/3D0 , 

1  ATHIRD=ASIXTH*2D0) 

S  =  EL 

IF  (EL.EQ.ODO)  GO  TO  2 
Cl  STARTER  BASED  ON  LAGRANGE'S  THEOREM 

G  =  IDO  -  Gl 
CL  =  DSQRT (IDO  +  EL**2) 

AL  =  DASINH (EL) 

W  =  G**2*AL/CL**3 
S  =  IDO  -  G/CL 

S  =  EL  +  G*AL/DCUBRT(S**3  +  W*EL*(1.5D0  -  G/0.75D0)) 

C2  TWO  ITERATIONS  (AT  MOST)  OF  HALLEY-THEN-NEWTON  PROCESS 

DO  1  ITER=1 , 2 

50  =  S*S 

51  =  SO  +  IDO 

52  =  DSQRT (SI) 

53  =  S1*S2 
FDD  =  G*S/S3 

FDDD  =  G* ( IDO  -  2D0*S0)/(S1*S3) 

IF  (ASIXTH*S0  +  Gl  .GE.  SW)  THEN 
F  -  (S  -  G* DASINH (S ) )  -  EL 
FD  =  IDO  -  G/S2 
ELSE 

F  =  SHMKEP (Gl ,  S)  -  EL 
FD  =  (SO/ (S2  +  IDO)  +  Gl)/S2 
END  IF 

DS  =  F*FD/ (AHALF*F*FDD  -  FD*FD) 

STEMP  =  S  +  DS 
IF  (STEMP. EQ.S)  GO  TO  2 

F  =  F  +  DS* (FD  +  AHALF*DS* (FDD  +  ATHIRD*DS*FDDD) ) 

FD  =  FD  +  DS* ( FDD  +  AHALF*DS*FDDD) 

1  S  =  STEMP  -  F/FD 

2  SHKEPL  =  S 
RETURN 
END 


Appendix  B 

AN  UNSOPHISTICATED  SHMKEP  PROCEDURE 

DOUBLE  PRECISION  FUNCTION  SHMKEP  (GI,  S) 

ACCURATE  COMPUTATION  OF  S  -  (1  -  Gl ) *DASINH (S ) 
WHEN  (Gl,  S)  IS  CLOSE  TO  (0,  0) 

IMPLICIT  DOUBLE  PRECISION  (A-H,  O-Z) 

G  =  IDO  -  Gl 

T  =  S/(1D0  +  DSQRT ( IDO  +  S*S) ) 

TSQ  =  T*T 
X  =  S* (Gl  +  G*TSQ) 

TERM  =  2D0*G*T 
TWOI1  =  IDO 
TWOI1  =  TWO I 1  +  2 DO 
TERM  =  TERM*TSQ 
XO  =  X 

X  =  X  -  TERM/ TWO I 1 

IF  (X.NE.XO)  GO  TO  1 

SHMKEP  =  X 

RETURN 

END 


Appendix  C 
THE  EKEPL  PROCEDURE 


DOUBLE  PRECISION  FUNCTION  EKEPL (EM,  El) 

C  KEPLER'S  EQUATION,  EM  =  EKEPL  -  (1  -  El) *DSIN ( EKEPL) , 

C  WITH  El  IN  RANGE  1  TO  0  INCLUSIVE,  SOLVED  ACCURATELY 

C  (BASED  ON  EKEPL3 ,  BUT  ENTERING  El  NOT  E) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

PARAMETER  (PI*3 . 14 1592653 5897932 384 6264 3 38 3 28 DO , TWOPI=2DO*PI , 

1  PINEG=-PI,  SW-9.25D0,  AHALF=0 . 5D0 ,  ATHIRD=AHALF/1 . 5D0) 

Cl  RANGE-REDUCE  EM  TO  LIE  IN  RANGE  -PI  TO  PI 

EMR  =  DMOD (EM, TWOPI ) 

IF  (EMR. LT. PINEG)  EMR  =  EMR  +  TWOPI 
IF  (EMR.GT.PI)  EMR  =  EMR  -  TWOPI 
EE  =  EMR 
IF  (EE)  1,4,2 

1  EE  =  -EE 

C  (EMR  IS  RANGE-REDUCED  EM  &  EE  IS  ABSOLUTE  VALUE  OF  EMR) 

C2  STARTER  BY  FIRST  SOLVING  CUBIC  EQUATION 

2  E  =  IDO  -  El 

W  =  DCBSOL (E ,  2D0*E1 ,  3D0*EE) 

C3  EFFECTIVELY  INTERPOLATE  IN  EMR  (ABSOLUTE  VALUE) 

EE  =  ( EE*EE  +  (PI  -  EE)*W)/PI 
IF  (EMR.LT.ODO)  EE  =  -EE 

C4  DO  TWO  ITERATIONS  OF  HALLEY,  EACH  FOLLOWED  BY  NEWTON 

DO  3  ITER=1 , 2 
FDD  =  E*DSIN (EE) 

FDDD  =  E*DCOS ( EE) 

IF  (EE*EE/6D0  +  El  .GE.  SW)  THEN 
F  =  (EE  -  FDD)  -  EMR 
FD  =  IDO  -  FDDD 
ELSE 

F  =  EMKEP ( El , EE)  -  EMR 
FD  =  2D0*E*DSIN (AHALF*EE) **2  +  El 
END  IF 

DEE  =  F*FD/ (AHALF*F*FDD  -  FD*FD) 

F  =  F  +  DEE* ( FD  *  AHALF*DEE* ( FDD  +  ATHIRD*DEE*FDDD) ) 

C*  TO  REDUCE  THE  DANGER  OF  UNDERFLOW  REPLACE  THE  LAST  LINE  BY 

C*  W  =  FD  +  AHALF*DEE* (FDD  +  ATHIRD*DEE*FDDD) 

FD  =  FD  +  DEE* (FDD  +  AHALF*DEE * FDDD) 

3  EE  =  EE  +  DEE  -  F/FD 

C*  IF  REPLACING  AS  ABOVE,  THEN  ALSO  REPLACE  THE  LAST  LINE  BY 

C*  3  EE  =  EE  -  (F  -  DEE* ( FD  -  W))/FD 
C5  RANGE-EXPAND 

4  EKEPL  =  EE  +  (EM  -  EMR) 

RETURN 

ENO 
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Fig  1  Maximum  relative  truncation  error  for  the  solution  procedure, 
SHKEPL,  for  the  reformulated  hyperbolic  equation  (e  <  10) 
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